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We carry out numerical simulations to investigate spontaneous vortex formation during a tem- 
perature quench of a superconductor film from the normal to the superconducting phase in the 
absence of an external magnetic field. Our results agree roughly with quantitative predictions of 
the flux trapping scenario: In fast quenches the agreement is almost perfect, but there appears to 
be some discrepancy in slower ones. In particular, our simulations demonstrate the crucial role the 
electromagnetic field plays in this phenomenon, making it very different from vortex formation in 
superfluids. Besides superconductor experiments, our findings also shed more light on the possible 
formation of cosmic defects in the early universe. 



I. INTRODUCTION 

When a superconductor in rapidly cooled from the 
normal to the superconducting phase in the absence of 
an external magnetic field, vortices (and antivortices) 
form. This phenomenon has been observed in several 



experiments 



1.2.3.4.5.6.7.8 



and has analogues in other con- 



densed matter systems such as superfluids 9 ! 10 ! 11 and 
liquid crystalsii2ii&i!ii£ In the latter cases, the or- 
der parameter is electrically neutral and defect forma- 
tion can be understood in terms of the Kibble-Zurck 
mechanism^£iiLi& In contrast, the Cooper pairs in su- 
perconductors are electrically charged, which means that 
the electromagnetic field plays an important role>i2i22i 

Understanding the mechanisms in which vortices form 
would also have implications well beyond condensed mat- 
ter physics. In particular, similar mechanisms may 
have produced topological defects such as cosmic strings, 
magnetic monopoles and domain walls in the early 
universe jiSiSi either in cosmological phase transitions 17 
or in brane collisions. 22 This is particularly interesting, 
because two possible gravitational lensing events by cos- 
mic strings have recently been observed ^S^k 

In cosmology, a neutral order parameter corresponds 
to a global symmetry and a charged order parameter 
to a local gauge symmetry. In the latter case, which 
is more common, the vortices are Nielsen-Olesen cosmic 
strings^ The role of the Cooper pairs is played by the 
Higgs field, which obtains a non-zero vacuum expectation 
value at the transition. 

In this paper, we investigate vortex formation in the 
case of a two-dimensional film in three-dimensional space. 
This is the relevant setup for actual superconductor ex- 
periments, but to our knowledge no simulations have 
been carried out in the past. While in the case of a neu- 
tral order parameter, a two-dimensional film can be stud- 
ied using a fully two-dimensional simulation, the three- 
dimensional nature of the electromagnetic field means 
that the same is not true for superconductors^ 7 - 

Because our setup captures some of the pecularities 



of fully three-dimensional systems, it is also more rel- 
evant for cosmological applications than previous two- 
dimensional studies. 28 A different borderline case be- 
tween two and three spatial dimensions has been studied 
in the past in Refs. HiJ andE^. 



II. SETUP 

The system we shall study consists of a two- 
dimensional superconductor film in a three-dimensional 
space. The superconductor is described by the Ginzburg- 
Landau theory with an electrically charged order param- 
eter that is confined to the film. In contrast, the electro- 
magnetic field lives in the full three-dimensional space. 
We shall denote three- vectors by bold-faced letters (such 
as x) and two-vectors on the film by letters with arrows 
on top (x). 

We consider a hypothetical experiment in which the 
whole system is first heated up to a temperature above 
T c> so that the film is in the normal phase. Then the 
system is gradually cooled through the phase transition 
to the superconducting phase. 

As usual, the two-dimensional Ginzburg-Landau (free) 
energy is written as 



E,i, = 



d 2 x 



I 

2m 



(-ihV - eAj 



The field ip(x,y) is a complex scalar field, and the two- 
dimensional vector potential A(x,y) is obtained from 
the full three-dimensional A(x, y, z) by restricting to the 
plane z = 0, and x and y directions, 

A x (x,y) = A x (x,y,0), A y (x, y) = A y (x, y, 0). (2) 

The energy of a three-dimensional electromagnetic field 
configuration is given by 



d A x 



e E 2 



2flQ 



(3) 
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where B = V x A is the magnetic field (or induction) and 
E is the electric field. If we choose the temporal gauge, 
we can write E = —dA/dt. In what follows, we shall 
simplify the notation by choosing a unit system in which 
m = 1/2, c = h = jjiQ = eo = Ub = 1- 

Because we are interested in the time evolution of the 
system, we need to include the kinetic energy term cor- 
responding to the time derivative n = ip, for which we 
choose the standard "relativistic" form. The total energy 
is then given by 



= / d 2 x 



M 2 + "H 2 + ^I 4 + I^I 2 



-~ / d 3 x(E 2 



B 



(4) 



one can say that it cools down. If a is negative, the 
system will eventually reach the critical temperature at 
which the film becomes superconducting. 

Note that Eq. J7J), with an appropriate rescaling, be- 
comes the usual time-dependent Ginzburg-Landau equa- 
tion in the limit a — > oo. However, a determines the rate 
of cooling of the system, which we want to be able to 
vary. Therefore, we have to keep the second-order terms 
both in Eq. Q and in the Maxwell equation. 

A damping term like this is not the only possible way 
to induce a phase transition. In Ref. Il9t this was done 
by varying the quadratic term a at constant tempera- 
ture. We do not believe the qualitative phenomena we 
are interested in are sensitive to the specific choice, but 
quantitative details may well be different. 



where we have used the covariant derivative Dip = (V 
ieA)ip. The full partition function is 



VipVirVAVE e- Etot/T . 



(5) 



This partition function describes the thermal equilibrium 
state of the system. We shall use the thermal ensemble as 
the ensemble of initial conditions for the time evolution, 
choosing a high enough temperature so that the film in 
initially in the normal phase. 

The equations of motion are obtained by using E tot 
as the Hamiltonian and identifying it* as the canonical 
momentum conjugate to ip. This gives the equations 



D 2 i> 



■ aip + f3\i>\ 2 ip 



0. 



A + VxVxA = j, 



(6) 



where j is the electric current. The latter equation is, of 
course, equivalent to Maxwell's equations. The^ current 
is confined on the plane z = 0, i.e., j = S(z)(j x , j y , 0), 

where the two-dimensional current is j = 2elnx0* Dip . 

Because Eq. (0 describes a thermal equilibrium state, 
the system stays in equilibrium if it is evolved according 
to Eqs. jnj. In order to cool the system, we modify the 
equations by adding a damping term to the equation for 



ip - D 2 ip + aip + (3\ip\ 2 ip = — crip. 



(7) 



In order not to violate Gauss's law V • E = p, we need to 
modify the equation for the vector potential A as well, 
by adding an Ohmic contribution to the electric current, 



j = S(z)0 x J v ,O) + aE, 



(8) 



where the damping rate a plays the role of the conduc- 
tivity. 

When the system is evolved according to these equa- 
tions, the total energy gradually decreases. Strictly 
speaking, the system will not be in thermal equilibrium, 
but using a suitable effective definition of temperature, 



III. FLUX TRAPPING 

As long as the system is in the normal phase, it will 
stay close to thermal equilibrium, because the film does 
not affect the dynamics of the vector potential signifi- 
cantly. To a good approcimation, the Fourier modes of 
the orthogonal magnetic field B = d x A y — d y A x with 
wave number k greater than a will oscillate with a de- 
creasing amplitude 



B(k) - cxp(-crt/2) cos(fct). 



(9) 



The equilibrium distribution given by Eq. (|SJ) is approx- 
imately Gaussian in the normal phase, and it can there- 
fore be completely characterized by the two-point func- 
tion 



(B(k)B(k')) = G(\k\)(2Tr) 2 S(k + k'), 



when 



,27 



G(k) 



Tk 
~2~' 



(10) 



(11) 



Thus, the modes with |fc| <C a retain the equilibrium dis- 
tribution with an exponentially decreasing effective tem- 
perature T ff (t) w Ti„i exp(— at). 

When the temperature reaches the critical value T c , the 
order parameter tp becomes non-zero and the dynamics 
becomes non-linear. If the system were to stay in equi- 
librium, it would now be in the superconducting phase 
and therefore repel magnetic fields. 

The equilibrium two-point function is^i 



G(k) 



Tk 



1 



2 1 + fcA' 



(12) 



where A is a screening length that starts from zero at T c 
and grows exponentially as the temperature decreases. 
At any finite temperature, thermal fluctuations with 
wavelength less than A are suppressed relative to the nor- 
mal phase, but longer wavelengths are unaffected. This 
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means that there is no sharp transition but the apparent 
critical temperature depends on the length scale. 

Even though it is impossible to solve the non-linear 
equations of motion analytically, we can generally say 
that the amplitude of a Fourier mode of B with wave 
number k cannot decay arbitrarily fast. Moreover, the 
longer wave lengths react slower. For instance, if the 
dynamics of the long-wavelength fluctuations is diffusive, 
the fastest possible decay rate is 7 ma x(&) = k 2 /D, where 
D is the diffusion constant. It is therefore unavoidable 
that if the transition takes place in a finite time, the 
longest wavelengths are too slow to react, and freeze out. 
This means that there is a critical wave number k c so 
that modes with k < k c still have approximately their 
initial amplitude after the transition. 

The survival of the long- wavelength fluctuations means 
that at distances less than X/k c , there is effectively a uni- 
form magnetic field. If this length scale is longer than the 
Pearl length, i.e., the size of a vortex, this magnetic field 
must form an Abrikosov vortex lattice M> This mechanism 
of vortex formation is called flux trapping. 

The number density of vortices per unit area produced 
by flux trapping is approximately^! 



e 

27 



2tt ' 



(13) 



Moreover, because the frozen-out magnetic field consists 
of modes with wavelengths longer than 2tt /k c , there are 
clusters of size 27r/fc c of equal-sign vortices. We can esti- 
mate the number of vortices in each cluster by assuming 
that they are disks of radius l/k c . 



U 2 T C 
87rfc c 



(14) 



vortices. For the estimate in Eq. 1(13(1 to be valid, N c \ 
has to be greater than one. Otherwise, one will have to 
consider the coupled dynamics of both the electromag- 
netic field and the order parameter. It seems likely that 
vortex formation can then be described by the Kibble- 
Zurek mechanism ,i£iiLi2i although the electromagnetic 
field may still play a role. 

In fact, the friction term a causes modes with \k\ < a 
to freeze out even in the absence of any critical dynamics. 
The solution of the linearized equations of motion with 
thermal initial conditions gives 



G(k) 



Te 



T 



"271? 

dk z P 
~2~7 k 2 



<J 2 , . a . , . 4k 2 

— cosh at + — sinn at — 

a z a a 1 



Tk 



e -(2k /a)t = v r { Cy /2kH/a, (15) 



where k 2 = k 2 + k 2 and a 2 = a 2 — 4k 2 , and Erfc is the 
complementary error integral. The second line is valid at 
late times when \k\ <C a. This gives rise to an apparent 
critical wave number k c — \J a jit. 



We are mainly interested in the physically more rele- 
vant case in which the freeze-out is caused by the critical 
dynamics, which restricts us to relatively low damping 
rates a. Nevertheless, simulations with higher a are also 
useful, because we can make more precise theoretical pre- 
dictions using the exact two-point function 1)15(1. 



IV. FLUX QUANTIZATION 

Let us now assume that at the freeze-out, the two-point 
function is given by some function G(k), and calculate 
how many vortices should be formed. If the amplitude 
of the frozen-out fluctuations is high enough, this can be 
done by considering a circular region of radius R. We 
assume that the magnetic field is more or less uniform at 
this length scale and that the magnetic flux $(i?) through 
this region is much greater than one flux quantum 27r/e. 
In that case, the number of vortices N(R) is given by the 
flux 



N{R)=nvR 2 = — \®{R)\ = — 



2tt 1 



2vr 



d 2 xB(x) 



, (16) 



where n is the number density of vortices per unit area. 
We can therefore write 



lim 



(t>(R) 2 ). 



(17) 



It is straightforward to write ($(i?) 2 } in terms of G(k), 



d 2 xd 2 y(B(x)B(y)) 



R 4 

R*_ 
4 



d 2 kG(k), 



(18) 



where J\{kR) is a Bessel function. Thus, we expect 

1/2 



e 
I^ 2 



d 2 kG{k) 



(19) 



In the friction-dominated case, we can combine the 
integrations in Eqs. l|15f) and (|19fl into one integral, 



e 2 Te 



12tt 4 
e 2 T 



rfkk 2 



a 2 , „ <r . , . 4k 2 

-^j cosh at + — sinn at ^5- 

a a a 



12tH 



dkk 2 e- 2k = 



e 2 T 



4S;r 2 Gvrt) 



3/2 



(20) 



where the approximate equality is valid at late times. In 
principle, t should be chosen to be the time the vortices 
form. 

For slower a, the two-point function G(k) cannot be 
calculated analytically, but we assume that it can still be 
approximated by a Gaussian function, 



G{k) 



2 



(21) 
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where T2 and fc 2 are constants. Then, Eq. (|19(l yields 

(22) 



e 2 T 2 fc 2 3 
64tt 5 / 2 ' 



The number of vortices per cluster is roughly given by 
iV c i w rnr/k 2 , which gives 



3 2 T, 



64 7 r 1 /2fc 2 ' 



(23) 



and in principle this number should be much greater than 
one for Eq. (|16|l . and consequently also these estimates, 
to be valid. 

It is important to note that Eqs. I|15l) and (|21|l are ide- 
alized descriptions of what the two-point function G{k) 
should look like at the time of the freeze-out, when the 
magnetic field is still smooth over long distances. When 
the vortices form, the magnetic flux is quantized and this 
introduces microscopic structure to the two-point func- 
tion. This will be important when we measure it and 
attempt to extract the fit parameters. 

Let us first assume that N vortices and N antivortices 
are formed, and that they are point-like so that their 
magnetic field can be described by a sum of delta fun- 
tions, 



2tt 



n 



(24) 



i=i 



where xf are the positions of the vortices and antivor- 
tices. The two-point function is 



Air 

B{x)B{y)) = —Y,G l] {x-y) 1 



(25) 



where the two- vortex correlator Gij (x — y) is 

G^x-y) = {[S(x-4)-5(x-x7)][6(y-x+)-6(y~x-)]), 

.(26) 

and the brackets () indicate integration over the positions 

As long as i ^ j, Gij{x — y) = Go(x — y) is independent 
of i and j. On the other hand, for i = j, we have 

G u (x-y) = (S(x-xf)S(y-xf)) 
+(5(x-xX)6{y-x7)) 
-{5(x-xf)5(y-x7)) 
-(S(x-xr)S(y-xt)). (27) 

The two last terms give contributions of order N/A 2 , 
which we will ignore but the two first give delta functions, 
so 



G u (x-y) = 25(x-y)/A, 
where A is the area of the film. 



(28) 



electric charge 


e 


0.3 


coupling constants 


a 


-0.25 




P 


0.18 


lattice size 


N 3 


512 3 


lattice spacing 


Sx 


1.0 


time step 


St 


0.05 


initial temperature 


Ti„i 


10.0 


thermalization cycles 


N th 


16 


time evolved in each 


tth 


32.0 



TABLE I: Parameter values 

In total, we have in Fourier space, 
4tt 2 



G(k) 



[(2N/A) + N(N - l)G {k)] 



(29) 



The continuous-field limit is basically equivalent to tak- 
ing TV — y 00, and 2N/A = n is the number density of 
vortices, so we find 



G(k) 



47T 2 71 



+ Gcont(k), 



(30) 



where G con t(fc) is the two-point function given by the 
continuous magnetic field. 

Of course, the magnetic field is not completely local- 
ized, which means that the delta functions spread over a 
finite distance. If we write the Fourier transform of the 
magnetic field profile of a vortex as Gi(fc), we have 

Air 2 n 

G{k) = — =-Gi(fc) + G cont (k). (31) 

We use the Gaussian ansatz in Eq. 121|) for G con t(fc)- 
In principle, G\(k) could be calculated by finding the 
static vortex solution, but we simply parameterize it by 
an exponential, so that the whole two-point function is 



G(k) 



2 



T2k (k/k^ 
2 



(32) 



We will use this form for fitting our numerical results. 
The above discussion predicts that k± should not depend 
on the nature of the quench, because it is a property of 
the static solution. The value of Xi should be propor- 
tional to the number of vortices. 

The linearized result in Eq. H15fl is not of exactly the 
same form as Eq. I|21[l. and therefore the parameters T2 
and fc 2 will not be exactly the same as their theoretical 
values in the fast-quench limit, but they should still be 
close to them. 



V. SIMULATIONS 

We tested our predictions by solving the fully non- 
linear equations of motion l|A5|l numerically using the 
initial conditions given by the partition function (jSJ). To 
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FIG. 1: Time evolution of the two-point function at a = 1.0. 
From top to bottom, the curves correspond to t = 0, 2, 4, 6, 
8, 10 and 20. 
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FIG. 2: The measured two point function G(k), and the 
linearized prediction calculated numerically from Eq. H5H 
[dashed lines]. From left to right, the curve pairs correspond 
to o = 0.25, a = 1.0 and a = 5.0. 



do this, we defined the system on a three-dimensional 
lattice, one slice of which corresponds to the film. The 
details of the lattice implementation are presented in Ap- 
pendix [S] 

In the preparation of the thermal initial conditions 
we employed the fact that the energy functional -Etot in 
Eq. (IA7D is quadratic and diagonal in E and it. There- 
fore, their probability distribution is Gaussian. The only 
complication is the Gauss law (|A6|) . without which the 
field values at different points would be uncorrelated. 

To generate the initial conditions, we first drew ran- 
dom, uncorrelated values for the component of tt that is 
parallel to tp, i.e., tt^ = ■0(Re , 0*7r)/(^*^). This compo- 
nent does not appear in the Gauss law (|A6|) . and there- 
fore it has the simple Gaussian probability distribution 



P ex P 



5x 2 2 



(33) 



Then, we went through all plaquettes and considered a 
change of the electric field at all links around the plaque- 
tte by the same amount, 



E 



E 

E ; 



i,(x) 



E i,(x) + e, 
Ej\(x)+j + e, 



E 



'j,(x) 



E, 
E 



(x)+i 

'i>(*0 ~ 



(34) 



Such a change does not change the divergence of E and 
therefore does not affect the Gauss law constraint. It 
would change the energy by an amount that is at most 
quadratic in e, 



A-Etot(e) = Qe 2 + Le, 



(35) 



where the constants Q and L depend on the field val- 
ues at the neighbouring links. Therefore the probability 
distribution for e is Gaussian, 



p [e] oc exp 



A^tot(e) 



T 



(36) 



It is straightforward to draw the value of e from this 
distribution. 

Next we evolved the field configuration for time tth us- 
ing the equations of motion (|A5|I , with zero conductivity. 
We repeated this cycle of randomization and evolution 
steps Nth times, monitoring the evolution of the two- 
point function G{k). When it reached the equilibrium 
form given by the analytical calculation, we considered 
the system to have thermalized. 

We then solved the equations of motion <| A5|> with the 
initial conditions produced by the thcrmalization algo- 
rithm. We used a number of different values for the con- 
ductivity a to test the dependence on the cooling rate, 
and for each value we repeated the run several times us- 
ing different initial conditions from the same ensemble. 
The values of the other parameters are shown in Table [I] 



VI. RESULTS 

Our main aim was to test the flux trapping theory, 
and therefore we measured the two-point function G{k) 
of the magnetic field fluctuations on the film. In Fig. ^ 
we show the time evolution of G(k) for a = 1.0. One can 
see that, as expected, the long- wavelength modes freeze 
out to an high amplitude, whereas short wavelengths are 
exponentially suppressed. The plot also shows how the 
amplitude at around 0.3 < k < 0.8 increases after t « 10, 
when the system enters the superconducting phase and 
the two-point function changes from Eq. (|21|) to Eq. I|32|) 
because of flux quantization. 
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FIG. 3: The exponential cutoff scale fci as a function of a. For 
o — > it approaches a constant value fci = 0.104(1), which 
characterizes the size of a static isolated vortex solution. At 
high a, vortex density is high and the system still away from 
equilibrium, and the vortices are not well described by the 
static solution. In this regime, we fit ki = 0.0140(6)o- l fi0(3) . 
In all other fits, we fixed fci = 0.1. 




FIG. 4: Two point functions G(k) together with fits of the 
form I32H . From top to bottom, the curves correspond to 
a = 5.0, 0.5 and 0.05. 



To make it possible to compare the results for dif- 
ferent values of cr, we chose to carry out our measure- 
ments at time t — 20/ cr, so that the effective temper- 
ature is the same in each case. Fig. [21 shows the mea- 
sured two-points function at that time for selected val- 
ues of a. They are compared with the results for lin- 
earized friction-dominated freeze-out, calculated numer- 
ically from the first line of Eq. (|15|) . The plot shows that 
the linearized result works very well in fast quenches. For 
a < 1, one can see a clear discrepancy, which is a sign 
that the non- linear effects have become important. The 
difference is also partly due to the contribution Gi(fc) 
from flux quantization. 

We fitted the two-point functions with the ansatz in 



Eq. (|32|l . In Fig. |3J we show k\ as a function of a. As 
expected, its value is independent of a in slow quenches. 
At higher cr, the vortex density becomes higher and non- 
equilibrium effects more important, and therefore it is not 
surprising that k\ starts to grow. This can be interpreted 
as the size of the vortex getting smaller as more of them 
are packed in the same area. 

To improve the accuracy of the fit parameters in sub- 
sequent fits, we fixed fci to 0.1. Some examples of these 
fits are shown in Fig. 0] The plateau at very short wave- 
lengths, which corresponds to modes still in thermal equi- 
librium was excluded from the fits. 

In Fig. 0we show the fit parameters T\, T2 and k% 
as functions of cr. The errors were estimated using the 
bootstrap method, and contain only the statistical error. 
We did not attempt to estimate the systematic error in 
any measurement due to the choice of the fitting function. 

In the high-cr limit, the two-point function should be 
given by Eq. I|15l) . and to the extent that it can be ap- 
proximated by a Gaussian, we expect Ti w T = 10 and 
&2 ~ \J a jit. This is confirmed by the measurements. 
In this limit the values of T\ show significant scattering, 
but this understandable, because G{k) is dominated by 
the Gaussian term and we had also fixed k\ to a very 
different value from the best fit. 

In the opposite limit of low cooling rates, the parame- 
ters seem to be well described by power laws, 



1.218(6) 



Ti - 9.76(9)<7 
T 2 ~ 2.6(5)cr a43 ( 6 ) 



0.186(5)a a318 ( 8 ). 



(37) 



The deviation from the linear prediction ki « a /2t is a 
sign of non-linear dynamics. It is interesting to note, but 
possibly a coincidence, that the behaviour of T2 appears 
to agree with the same power law even at a ~ 10. 

Fig. [S] shows the number of vortices plus antivortices 
measured at time at = 20. The agreement with the an- 
alytical prediction Eq. (|2U|) shown by the dashed line is 
good for fast quenches cr > 1, apart from a constant 
factor of about 2. This indicates that vortices are pre- 
dominantly formed by flux trapping. 

For a < 1, we can fit the data very well with a power 
law N = 772(18)a 1 - 20(2) . This is compatible with the 
expectation that T\ oc N. We have also plotted the the- 
oretical prediction in Eq. (|22|l calculated using the fitted 
parameter values. As the plot shows, it gives the correct 
order of magnitude, although it does not agree perfectly 
for slow quenches. The prediction seems to suggest a 
different power-law behaviour from the observed one. 

As Fig. shows, the fluctuations have frozen out, but 
it is possible that their amplitude is not high enough to 
actually dominate the process of vortex formation. In- 
deed, if one estimates the typical number of vortices in 
a cluster using Eq. I|23[) one finds it decreases rapidly at 
low a (see Fig. 0. In principle, the flux trapping mech- 
anism requires this number to be greater than one. This 



7 




FIG. 5: Fit parameters T%, T2 and &2 as functions of a. (a) Coefficients T\ (filled circles) and Ti (empty diamonds) of 
the exponential and Gaussian terms in Eq. 13211 respectively. The dotted line shows the power-law fit T\ — 9.76(9)cr 1,218 ' 6 - 1 
to the first six data points. The dashed line shows a fit Tn = 2.6(5)<t°' 43 ' 6 - ) to the first four data points, (b) The critical 
wave- number k 2 . The solid line shows the theoretical predition k 2 = ^fa/2t, and the dashed and dotted lines are power- law 
fits k 2 = 0.186(5)<r°' 318(8) and k 2 = 0.115(7)<t°' 95(3) of the first five and last six data points, respectively. 




FIG. 6: Vortex number measured at at — 20 as a function 
of a. The dashed line shows the theoretical prediction I20H 
for fast quenches, and the solid line a a power-law fit N = 
772(18)<7 1 ' 20(2) of data at a < 0.5. The crosses show the 
prediction 12211 for slow quenches based on the fit parameters 
T 2 and k 2 . 



FIG. 7: The number of vortices in each cluster calculated 
from the measured vortex density and the fit parameter k 2 . 



typical number of vortices in a cluster. 



suggest that the slow quenches could possibly be better 
described by the Kibble-Zurek mechanism. 

It is possible to increase N c \ by using different pa- 
rameter values, and we carried out one set of runs with 
T = 1000. The clustering is indeed evident in the snap- 
shot of the magnetic field at time t — 400 in a run with 
a = 0.25 shown in Fig. [S] In the same figure, we have 
also plotted the quantity, nc(r) which measures the aver- 
age winding number around a circle of radius r centered 
at a vortex^S, Clustering is indicated by values greater 
than one, and the maximum value gives a measure of the 



VII. CONCLUSIONS 

In fast quenches a 3> 1, the two-point function G(k) 
of the magnetic field fluctuations is consistent with the 
predictions of the linearized theory. This allows us to 
calculate the number of vortices produced by flux trap- 
ping, and this prediction agrees very well with the mea- 
sured values. This provides strong support for the flux 
trapping scenario of vortex formation. Furthermore, the 
spatial distribution of vortices shows the tell-tale signal 
of vortex clusters. 
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APPENDIX A: LATTICE DISCRETIZATION 



FIG. 8: The magnetic field configuration on the film at time 
t = 400 in a run with T = 1000 and a = 0.25. The curve in the 
inset shows the quantity nc{r), which measures clustering. 



It should be kept in mind that in fast quenches, the 
relevant degrees of freedom are overdamped. This is not 
necessarily a problem, because it can be given a physi- 
cal interpretation as a relatively high conductivity. Fur- 
thermore, overdamped dynamics is commonly used in 
studies of superconductors and vortex formation^ Nev- 
ertheless, for many applications underdamped dynamics 
would probably be better. Ideally one should carry out 
the simulation without introducing damping at all^Si 
but even in our case, the relevant degrees of freedom are 
underdamped in slow quenches. 

The underdamped nature of dynamics makes it much 
more difficult to describe slow quenches theoretically, and 
therefore the predictions are not as robust. In spite of 
this, our results show a reasonably good agreement be- 
tween theory and simulation, although the scaling of the 
vortex number with the cooling rate a does not seem to 
agree. This may be an indication that the amplitude of 
trapped fluctuations is so low that vortex formation can 
be better described by the Kibble-Zurek mechanism. 

Simulations with a higher initial temperature support 
this picture, and show clear vortex clustering even in slow 
quenches. This is an unmistakable sign of flux trapping 

To understand better how vortices form when the am- 
plitude is not high enough for flux trapping, it would be 
important to derive theoretical predictions of the Kibble- 
Zurek scenario and compare them with our results. In 
particular, one must understand how a weak background 
magnetic field biases the Kibble-Zurek mechanism. 

Our results show that under the right conditions, the 
electromagnetic field plays an important role in defect 
formation. It remains to be seen if these conditions are 
satisfied in actual superconductor experiments. Likewise, 
the cosmological consequences of this remain largely un- 
explored. 



To write down the discretized equation of motion, we 
define the forward and backward derivatives 



Af / (X ) = ±Sx~ 1 (/ (x±J) - / (x) ) 



(Al) 



where Sx is the lattice spacing and i is a unit vector in 
the i direction. Similarly, we define the time derivative 

At/ (t ) = Sr 1 (/ (t) - f (t _ st} ) , (A2) 

where St is the time step. Using the link variable 

Ui = exp (ieSxAA , (A3) 

we also define the corresponding covariant derivatives on 
the film 



~1p(x) = Sx 1 (tfi^3)i>(g+t) - -0(5)) , 
~ip {s) = Sx- 1 (y {3) - Uf^_ t) ip^_ t ^ . (A4) 

The discretized equations of motion are 



DP 



AtEj^x) = e ijk£klmAj A+ A mi ( t)X ) - &Ei^t,x) 



jklm 

2e 



5 -8 z l^l t3) D l V(t,2), 

A t 7T( t>i ?) = ^ D~D?lf>( t ,£) - (JTT(t.S) 

i 

-aip( t .s) ~ P\ip(t.x)\ 2 ^(t,x), 



(A5) 



where the summation over directions goes from 1 to ei- 
ther 2 or 3 depending on the context. 
The lattice version of the Gauss law is 



2e 



(t,x) 



—Szhmp^n^a). (A6) 



This equality is exactly conserved by the lattice equations 
of motion and is a constraint the initial field configuration 
must satisfy. 
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The lattice version of the energy functional Q is 



jk 



Y,Sx 2 \tt*tt - —J2Rei>* { g } U h(S) 4> {S+i) 



4 
Sx 2 



n f ! + • (ATI 



It should be noted that both the equations of motion 
(|A5I) and the energy functional l|A7[) are gauge invariant. 
In this non-compact lattice formulation, the gauge group 
is, in fact, R rather than U(l), which has the advantage 
that vortices are topologically stable even on lattice. 

In the actual simulations, we used a finite lattice with 
periodic boundary conditions in all three directions. This 
has the consequence that the net magnetic flux through 
the film vanishes. 
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